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Abstract 

A parabolic equation for the propagation of periodic internal waves over varying 
bottom topography is derived using the multiple-scale perturbation method. Some 
computational aspects of the numerical implementation are discussed. The results 
of numerical experiments on propagation of an incident plane wave over a circular- 
type shoal are presented in comparison with the analytical result, based on Born 
approximation . 
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1 Introduction 



The parabolic equation method, as an approximation to some elliptic prob- 
lems, has been extensively used in mathematical physics. It was introduced 
by Leontovich & Fock [T] in the theory of electromagnetic wave propagation, 
applied in monograph of Babich & Buldyrev [2] to various diffraction prob- 
lems and developed by Tappert [3] and his successors in underwater acoustics. 
Later it appeared also in the theory of surface water wave propagation in 
works of Liu & Mei [3], Radder [5], Kirby Sz Dalrymple [B] and others. Since 
then, it has been proven through the contribution of several authors to be an 
effective model for dealing rapidly and accurately with propagation problems 
in coastal areas. 
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The standard, or narrow angle parabolic wave equation has the form of the 
quantum mechanical non-stationary Schrodinger equation and describes the 
propagation of waves in weakly inhomogeneous media, at small angles with a 
preferred direction (taken to be in this paper the x-direction). Being of the 
evolution type equation, it can be solved in the frame of the initial-boundary 
value problems and has an advantage over the geometric optic method in 
describing the wave field near a caustic. In the theory of the surface water 
wave propagation the parabolic equation method now can be considered as 
classic and is well explained in Mei's monograph [7]. In this theory some more 
complicated models, which describe wide-angle propagation [S] and include 
some nonlinear effects [6119] . were also developed. 

In the theory of internal wave propagation, though the geometric optic method 
is well-developed [TU|11II12] , little is known about the parabolic equation method, 
except the Kadomtsev-Petviashvili equation, which in the variable topography 
case was derived for the interfacial waves in [T5] . 

The case of continuous stratification is in close analogy with the acoustic 
case, where the adiabatic mode parabolic equations were first derived by the 
factorization method in [T3j and later by the multiple-scale method in [15] and 



The aim of our work is to obtain a standard parabolic equation for the internal 
wave propagation in the most simple cases and briefly discuss related compu- 
tational aspects. Rotation, which is ignored at this stage, can be included in 
the model as weak rotation [13], and will be considered in our further publi- 
cations. 

As an illustration we present some numerical calculations performed in the 
case of refraction of an incident plane wave on a circular-like shoal. 



2 Formulation and scaling 

The system of linear equations, describing small amplitude motions of strat- 
ified inviscid fluid with harmonic dependence on time t by the factor e~ luJt , 
from which we are starting, is 



[16] 



— iupou + P x = 
-iujp Q v + Py = 



-iujpow + P z + gpx = 




-lupx + wp 0z = 
u x + v y + w z = 



2 



Here x, y, and z are the Cartesian coordinates (vertical axis is directed up- 
ward), p — Po( z ) is t ne undisturbed density, P is the pressure, g is the gravity 
acceleration, pi = pi(x, y, z) is the perturbation of density due to motion, and 
u, v and w are respectively the x, y and z components of velocity. 

We consider these equations with the boundary conditions 

w = at z = 
w + uH x + vH y = at z = -H 

where H = H(x, y) is the bottom depth. 

To begin the multiple-scale procedure [17], we introduce a small parameter e, 
the slow variables X = ex and Y = e 1 / 2 ?/, the fast variable £ = (l/e)Q(X,Y) 
and expand the dependent variables as follows: 

u = u + eui + . . . , v — e 1/2 f 1/2 + . . . , 

w = w + ew\ + • • • , P = Po + ePi + . . . 

The scaling used in the definition of the slow variables is characteristic for the 
parabolic equation method [TUTl Section 4.10], with x-direction as the principal 
propagation direction and y-direction as the transverse direction. 

From now on we assume that 

H = H (X) + eH 1 (X,Y), 

so H is independent of Y. 

We also expand p\ in powers of e 

pi = Pio + epn + • • • ■ 



Substitution of these expansions into system ([T]) and the boundary conditions 
§2§ and changing the partial derivatives for the prolonged ones by the rules 



leads to the system of equations 

d 1 d 

-iou Po (u + e Ul + ...) + e(— + -Q X —)(P + ePi + ...) = , (3) 

-ia;p (e 1/ ^i/ 2 + ...) + e 1/2 ( ^ + V4)(fl> + ePi + •••) = , (4) 

-icjpoOo + e^i + • • •) + + + • • • + #(pio + epn +...) = > ( 5 ) 

-iw(pio + epu + •••) + (™o + erni + ■ ■ -)Poz = , (6) 



e( ^ + 7 4 )(Mo + eMl + --° 

+ el/2 (l^ + -@y^)(e 1/2 ^i/2 + •••)+ wo* + ™u + • • • = , 
with the boundary conditions 

Wo + ewi + . . . = at z — , 

wo+ewi + . . . + e(w + euj + . . -)g^( H o + tH i) 

+ e 1 / 2 (e 1 / 2 v 1/2 + ...)-^(H + eH 1 )=0 at z = -(H + eH r ) . 



(7) 



(9) 



3 The parabolic equation 



Equating coefficients in Eqs. (EH) of like powers of e, we obtain equations 
describing u h w h v i+1 / 2 , 1 = 0,1,.... 

At the order —1/2 in e we obtain from (pO) 

e y p 0§ = o , 

and put 0y = 0, so in the sequel depends only on X, = Q(X). 

At the zeroth order in e we have 

- iup u + Q x Po^ = , (10) 

- iupowo + P 0z = -gpxo , (11) 
- iup w + w poz = , (12) 

Q X U0!;+W0z = 0, (13) 

with the boundary conditions 

wq = at z = , 
wq = at z = —H . 
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After substitution in Eq. (iTTj) pxo from Eq. ffT2l we get 

O 2 po + gp 0z )w + iujP 0z = . (14) 
Twice differentiating this equation with respect to £ and using the expression 

= \ 2 MA) w o*)«, 
\Px) 

obtained from Eqs. ffTUj) and ffT3l) . we get 

(6x) 2 (w 2 po + gpoz)w tf + UJ 2 (p0zW0z)z = . 

We seek a solution of this equation in the form of WKB-type anzatz = 
A(X, Y)<f)(X, z)e^, where <fi is an eigenfunction with the eigenvalue k 2 = (©x) 2 
of the spectral problem 



u 2 (po4> z )z - k 2 (u 2 p + gp 0z )4> = 
0(O)=0(-tf o ) = O, 



(15) 



normalized by the condition 



Po ■ (<Pz) 2 dz = - (uj 2 p + gpo z )4> 2 dz = 1 . (16) 

This problem is known as the main spectral problem of the linear internal 
wave theory. 

At 0(e 1 ^ 2 ) we have only one equation 

- iupoVx/2 + Poy = , (17) 

which express the balance in the transverse direction for the quantities of order 
< 0(e). 

The system of equations at the first order in e is 

- iup oUl + 6 X P 1? + Pqx = , (18) 

- iupoW! + P Xz = -gpn , (19) 

UOX + @X^1? + V1/2Y + w lz = , (20) 

- vjjpxx + wiPoz = , (21) 
with the boundary conditions 

Wl = at z = 0, (22) 

and 

Wq(z) + ewi(z) + eu (z)H x = at z — —H — eHi . 
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Expanding velocities in Taylor series with respect to z at z = —Ho and col- 
lecting terms at e 1 , we reduce the last boundary condition on z = —H 

w\ - wq z Hi + uqHqx = at z = -H . (23) 



Considerations, similar to those in the derivation of the spectral problem 
Eq. ([T5]) (the details are given in Appendix |AJ, lead at 0(e) to the equation 

U) 2 

(u 2 po + gp z)w 1 ^+-^(p w lz ) z 

( , ^ 1 r, kx, A (24) 

Seeking solutions of Eq. ff24|) which depend on £ by the factor exp(i£), we 
obtain for w± the equation 

^ 2 (poWi 2 ) 2 - k 2 (u 2 p Q +gp Qz )wi 

= -uj 2 \ 2(pqUq X )z + T^--Poyy 2 - -r-(poUo)z J • ^ 
\ iu; /e / 

A solution to this differential equation with respect to z with the boundary 
conditions Eqs. (1221) and (1231) can be found only when a certain compatibility 
condition is satisfied, because the right hand side of Eqs. (125!) and ( 1231) contain 
the solution of the spectral problem (fl~5]) . The compatibility condition yields 
the required evolution equation for the amplitude function A 

A x + ^tA yy - \^A + ^Po(-Ho) ■ H ox ■ (MX, -H )) 2 A 

2ik 2 k 2k 2 (26) 

- 7^Po(-H ) • H x ■ (cf> z {X, -H )) 2 A = , 

which we call the parabolic equation for periodic internal waves. Rewritten in 
the initial coordinates [x,y), it has the form 

Ax + 2^ Ayy ~ + 2% Pq{ ~ Hq) ' H ° x ' ~ Ho))2 A 27 

- 7^Po(-Ho) ■ H x ■ (<p z (x, -H )) 2 A = , 

where Hi(x,y) = eiTi(ex, e 1 / 2 ?/), For detailed derivation of Eq. (126]) see Ap- 
pendix [HI 

For the numerical calculation of the coefficients of Eq. (126]) can be used, in 
principle, any algorithm for solving the spectral problem Eq. (fT5|) . Some prob- 
lems arise with the derivatives, namely, kx and 4> Z {X, —H j). 
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To avoid numerical differentiation in the calculation of kx we differentiate 
the spectral problem Eq. ([15]) with respect to X and get the boundary value 
problem for (f>x 

uJ 2 (p (j)xz)z ~ k 2 (u 2 p + gp 0z )(px ~ 2kk x (uj 2 p + gp 0z )^ = , 
(f>x(X, 0) = , <j> x (X, -H ) - H xMX, -H ) = . 

The compatibility condition for this problem is 



kx uj 2 



k 2k 2 



Po{-Ho)-Hox-{<p z {X,-Ho)y 



which gives the stable formula for the numerical calculation of kx/k (modulo 
the calculation of (f> z (X, —Ho))). 

Now we derive a formula for the stable calculation of the derivative <p z {X, —Ho). 
To do this, we multiply Eq. (fT5l) by z and integrate from — Ho to zero. After 
integrating by parts and some transformations we get the required formula 

1 / k 2 f° f° \ 

<f>z(X, -H ) = ~5 / ( uj2 Po + 9Poz)<f>- zdz- p 0z <f>dz ]. 

Po{-Ho)H \uj 2 J-Hq J-h J 



4 Numerical experiments 



The numerical experiments were conducted for a fluid with an exponential 
density stratification p = exp(— 72), where 7 > 0. In this case the complete 
set of normalized solutions of the spectral problem Eq. (fT5|) is 



CO 



2H \ 



7 2 # 2 + 4n% 2 



91 



10* 



4> n (z) = C exp (jz/2) sin (rnrz/H ) 
where C is the normalizing constant 



From Eqs. 



and we have also 

4>nz{-H ) = 



nn 



n 



n 



1)»C7— exp(- 7 #o/2). 
Ho 



(28) 



(29) 



For the numerical experiments the problem of the scattering of the incident 
plane wave of a given mode on the localized small inhomogeneities of the bot- 
tom topography was considered. So Ho is taken to be a constant and H\ (x, y) 
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Fig. 1. Bottom topography forms of the compact inhomogeneities (130p used for the 
numerical simulation. For all cases hi = lm, a = 500 m. (a) M = 0, Qo = 0; 
(b) M = 1, a = 0; (c) M = 1, a = vr/2; (d) M = 3, a = 0. 

describes the inhomogeneity, the total depth is Hq + Hi(x,y). For the inho- 
mogeneities of the form 



H x (r, a) = -A M r M cos(M(a + a )) exp(-r7(T 2 ) , M = 1, 



(30) 



where (r, a) are the polar coordinates centered at the point (x ,y = 0) with 
a = corresponding to the positive x-direction, the scattering problem for 
the propagating in x-direction incident plane wave of the nth mode with the 
vertical velocity w inc = A exp(ik n (x — x ))(p n admits an approximate solution 
in which the nth mode component of the scattered field is 



w s n cat (r, a, z) =A n (r } a) 



x exp( 



AnA M G i M+1 



lO^M' 

kt cos a) 



._ K M a 2M+2 e i(fc„r-vr/4) 



2 MH 



(31) 



) cos M(ip — aio) cos a) ■ 4> r 



where G = uj 2 p(—H 



—Hq)) , k = JZk 2 — 2k\ cos a and tan(^) 



sin a I (1 — cos a). The quantity that will be compared with the solution of the 
parabolic equation (1271) is the absolute value of the nth mode part amplitude 
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Fig. 2. Transverse cross-sections at x = 12000 m of relative wave field amplitudes 
in scattering of the 2nd mode internal waves on the shoal (|30|) with parameters 
M = 0, a = 500 m, «o = (Fig. QJa)). — , the parabolic equation method; — , the 
Born approximation (I3ip . 

of the incident + scattered field 

\A{x,y)\ = \A exp(ik n (x - x )) + A n (^J (x - x ) 2 + y 2 , arctan(y/ (x - x )))| . 



The solution Eq. ( 13TT) was obtained by the methods of the work [18] in the 
frame of the Born and far field approximations. The methods of the work [18] 
do not use any assumptions on the preferred propagation direction and are 
closely related to the methods of the work [T9] . 

The numerical experiments were conducted on the computational domain 
[0 < x < 12000] x [—5000 < y < 5000], where x and y are measured in 
meters. Eq. fl27j) with the initial condition A = A(0, y) = constant was inte- 
grated on the grid 395 x 349 by the Crank-Nicholson scheme [20J . In order to 
exclude reflections at the boundaries, the absorbing Baskakov-Popov bound- 
ary conditions were used |21j . which were adapted for non- vanishing at the 
boundaries initial conditions. 

The inhomogeneities in all cases have width parameter a = 500 m and the 
amplitude h\ = Am exp(— M/2) laJM/2 1 = 1 m with the center positioned 
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Fig. 3. Transverse cross-sections at x = 12000 m of relative wave field amplitudes 
in scattering of the 2nd mode internal waves on the shoal (|30p with parameters 
M = l,a = 500 m, «o = (Fig. U[b)). — , the parabolic equation method; — , the 
Born approximation (|3ip . 

at Xo = 2000 m, yo = m. The shape parameters M and ao were varied (see 
Fig. [p. The depth of the flat bottom, H , was taken to be 60 m. The density 
stratification parameter 7 = 0.00025, which corresponds to the Brunt- Vaisala 
frequency N = ^fg^ « 0.05 s~ x . 

The computations were done for the 1st mode and 2nd mode incident plane 
waves with 15 min time period. The transverse cross-sections of the 2nd mode 
computed field amplitude at x — 12000 m are presented in Fig. [2}l5]in compari- 
son with the Born type approximation scattering results obtained by Eq. ( I3TT) . 
The results for the 1st mode are analogous. 

Considering the results of computations, it is worth noting that the solution 
(13TT) has an approximative character, so the comparison with it is not exactly 
the test of accuracy of the derived parabolic equation. Nevertheless, since 
Eq. fE?T|) is of quite different genesis, and, in particular, free from any assump- 
tions on the preferred propagation direction, this comparison can lead to the 
conclusion that the parabolic equation describes sufficiently well the waves 
with propagation angles up to 45° (see the discussion on the propagation 
angles in [5]), scattering on the enough rough topography. 
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Fig. 4. Transverse cross-sections at x = 12000 m of relative wave field amplitudes 
in scattering of the 2nd mode internal waves on the shoal (|30p with parameters 
M = 1, a = 500 m, «o = tt/2 (Fig. He)). — , the parabolic equation method; — , 
the Born approximation (I3ip . 



5 Conclusion 



For the propagation of periodic internal waves over uneven bottom topography 
with small irregularities, the narrow-angle parabolic equation ( 1271) has been 
derived. It also takes into account slow, but not necessary small, variations of 
bottom topography in principal propagation direction. 



We have illustrated the use of the obtained equation by presenting the results 
of scattering of the plane internal waves over shoals of the special forms (I5U|) . 
The results of computations are in a sufficiently good agreement with the 
analytical solution (I3ip . obtained in the frame of the Born approximation, 
and support the applicability of equation (J27J) for computing of internal wave 
fields over uneven bottom with restrictions typical for the parabolic equation 
method in general [3lT7j . 



11 



1.05 




095 1 1 1 1 1 1 1 1 1 1 1 

-5000 -4000 -3000 -2000 -1000 1000 2000 3000 4000 5000 

y (m) 

Fig. 5. Transverse cross-sections at x = 12000 m of relative wave field amplitudes 
in scattering of the 2nd mode internal waves on the shoal (|30p with parameters 
M = 3, a = 500 m, «o = (Fig. [3Jd)). — , the parabolic equation method; — , the 
Born approximation (I3ip . 

A Derivation of Eq. (124ft 



From (jT9ll2Tj) we obtain 

(uj 2 p + gp 0z )w 1 + iuP lz = . (A.l) 

From ([18]) we get 

Pi& = ■^-iupoum - ■pz-Poxt ■ (A.2) 
From ([20]) . taking into account ([17]) . we have 

1 11 

«1£ = -T^— {Uq X + V 1/2 Y + % = -7^— [UOX + -. PoYY + W U ) ■ (A. 3) 

From (1101) we get 

P xz = (j^-PoUo^j = -iuj^-pou + iuj^-p u ox ■ (A.4) 
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Substitution of and Pqx^ from Eqs. (1A.31 IA.4[) into Eq. (1A.2j) gives 

-^iK = -7? iw Po ( M ox + r^— Pqyy + w lz ] 
k 2 \ iup J 

- -j^ (-iuj-^-p u + iujp u ox ) (A.5) 

If 1 k x \ 

= ~-gi luJ 1 2 p Mox + ^ p oyy + Po w iz - ~jrPo u o J ■ 

Differentiating Eq. (1A.5I) with respect to z and substituting the result into the 
second ^-derivative of Eq. (lA.ip . we obtain Eq. ff24]) . 



B Derivation of the compatibility condition 



In this appendix the compatibility condition for the boundary value problem 
Eqs. (E3 [22] and [23]) is derived. 



From Eq. f[T4"]) we have 



— Pqyyz = -^{^ 2 po + gPoz)w OY Y ■ 



IUJ 



UJ' 



From Eq. ffT3]) . taking into account that u and w depend on £ by the factor 
e 1 ^, we have 

1 i 

m ° = ~~k W ° z ' ° r U ° = k W ° z ' 
Substitution of these expressions into the right hand side of ( [25]) (denote it by 
RHS) yields 



RHS = - uj 2 



2j;{poWozx)z - 2^k x {pow 0z ) z + \ (uj 2 p + gp 

K K UJ 



z Wqyy 



ik x / 



-UJ 



2^A x (p (f) z ) z + 2^A(p (j)xz 



+^ (uj 2 Po + gpoz) A YY <\> - 3 1 -j^-A(p, 



Wz z 



■ exp(i£) 



fB.l} 



To obtain the required compatibility condition we multiply Eq. (1251) by the 
eigenfunction and integrate with respect to z from — H to 0. Then twice 
integrating by parts Eq. (125}) using Eqs. (122] [23]) we obtain 



f° RHS ■ (j)dz = uj 2 p (—H )a(j) z (X, —H ) 

J -H 



(B.2) 
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where 

a =Wx(X, -H ) = w 0z (X, -H ) ■ H x - u (X, -H ) ■ H ox 

= (H, - l -Hox)wo z {X, -Ho) = (H x - ~H 0X )A<i> z (X, -H ) ■ exp(i£) . 

From Eq. flRll and Eq. flB~2j) we get 



H 



{po4>z)z4>dz + \a yy I (u 2 p + gpoz)4> A dz 



Ho 



,(X,-H )) 2 . 



k x f° i f° 

-3i—A (p (p z ) z (j)dz + -A 2(p (fixz)z(f>dz 

rv J — Ho fv J—Hq 

= - P o{-Ho)-{H l - l -Hox)A- 
Using the normalizing condition Eq. (JT6l) we have 

{^j H0 P0{(t)z)2dZ ) x 

= -2-^ + ^ £ h 2po<p z <Pxz dz - po(-Ho) ■ (MX, -H )) 2 H ox ^ = 
and thus 



(B-3) 



k J-Ho 



2(po4> Xz )z4>dz 



k J-Ho 



2p 



0<PXz<Pz 



dz 



-yApo(-Ho) ■ (MX, -Ho)fHox - 2i^A. 



k 



With this and the normalizing condition Eq. f|T6|) we obtain from Eq. (IB.3I) 
the required parabolic equation 

1 1 / 2 

Ax + ^A YY - - -fA + ^Po(-Ho) ■ Hox ■ (MX, -Ho)) 2 A 

- ^Po(-Ho) ■ Hx ■ (MX, -Ho)f A = . 



References 



[1] Leontovich, M. A., Fock, V. A., 1965. Solution of the problem of propagation 
of electromagnetic waves along the Earth's surface by the method of parabolic 
equations. Chapter 11, Electromagnetic Diffraction and Propagation Problems, 
ed. V. A. Fock. Pergamon Press. 

[2] Babich, V. M., Buldyrev, V. S., 1972. Asymptotic methods in problems of short- 
wave diffraction. Nauka, Moscow. (In Russian). 



14 



[3] Tappert, F. D., 1977. The parabolic approximation method, in Wave 
Propagation and Underwater Acoustics, ed. J. B. Keller, J. S. Papadakis. 
Lecture Notes in Physics, Springer- Verlag, Berlin and New- York. 

[4] Liu, P. L. F., Mei, C. C, 1976. Water motion on a beach in the presence of a 
breakwater. 1. J. Geoph. Res. 81, 3079-3094. 

[5] Radder, A. C, 1979. On the parabolic equation method for water-wave 
propagation. J. Fluid Mech. 95, 159-176. 

[6] Kirby, J. T., Dalrymple, R. A., 1983. A parabolic equation for the combined 
refraction-diffraction of Stokes waves by mildly varying topography. J. Fluid 
Mech. 136, 453-466. 

[7] Mei, C. C, 1983, The applied dynamics of ocean surface waves. John Wiley & 
Sons, New York et al.. 

[8] Dalrymple, R. A, Kirby, J. T., 1988. Models for very wide-angle water waves 
and wave diffraction. J. Fluid Mech. 192, 33-50. 

[9] Liu, P. L.-F, Yoon, S. B, 1985. Nonlinear refraction-diffraction of waves in 
shallow water. J. Fluid Mech. 153, 185-201. 

[10] Keller, J. B., van Mow, C, 1969. Internal wave propagation in an inhomofeneous 
fluid of non-uniform depth. J. Fluid Mech. 38, 365-374. 

[11] Miropolsky, Yu. Z., 1974. Propagation of internal waves in theocean with 
horizontal inhomogeneities of the density field. Izvestiya USSR Academy of 
Sciences, Physics of Atmosphere and Oceans 10(5), 519-532 (in Russian). 

[12] Voronovich, A. G., 1976. The propagation of surface and internalwaves in an 
approach of geometrical optics. Izvestiya USSR Academy of Sciences, Physics 
of Atmosphere and Oceans 12 (8), 850-857, (in Russian). 

[13] Chen, Y., Liu, P. L.-F., The Kadomtsev-Petviashvili equation for interfacial 
waves, 1995. J. Fluid Mech. 288, 383-408. 

[14] Collins, M. D., 1993. The adiabatic mode parabolic equation. J. Acoust. Soc. 
Amer. 94, 2269-2278. 

[15] Trofimov, M. Yu., 1999. Narrow-Angle Parabolic Equations of Adiabatic Single- 
Mode Propagation in a Horizontally Inhomogeneous Shallow Sea. Acoust. Phys. 
45, 575-580. 

[16] Trofimov, M. Yu., 2002. Wide-Angle Adiabatic Mode Parabolic Equations. 
Acoust. Phys. 48, 728-734. 

[17] Nayfeh, A. H., 1973. Perturbation methods. John Wiley &: Sons, New York et 
al.. 

[18] Zakharenko, A. D., 2002. Scattering of internal waves from small sea bottom 
inhomogeneities. Proceedings of The Sixth Pan Ocean Remote Sensing 
Conference ( PORSEC), Bali, 3-6 September 2002, Vol. II, pp. 773-777, 
|arXiv:physics/070122Tj 



15 



[19] Llewellin Smith, S. G., 2002. Conversion of the barotropic tide. J. Phys. 
Oceanogr. 32, 1554-1566. 

[20] Potter, D., 1973. Computational physics. John Wiley & Sons, N.-Y. 

[21] Baskakov, V. A., Popov A. V., 1991. Implementation of transparent boundaries 
for numerical solution of the Schrodinger equation. Wave Motion 14, 123-128. 



16 



